AN ANALYSIS OF BROKEN Pi-NONCONFORMING FINITE ELEMENT 
METHOD FOR INTERFACE PROBLEMS 

DO Y. KWAK* AND K. T. WEE+ 

Abstract. We study some numerical methods for solving second order elliptic problem with interface. We introduce 
an immersed interface finite element method based on the 'broken' Pi-nonconforming piecewise linear polynomials on 
interface triangular elements having edge averages as degrees of freedom. This linear polynomials are broken to match 
the homogeneous jump condition along the interface which is allowed to cut through the element. We prove optimal 
orders of convergence in H 1 and L 2 -norm. Next we propose a mixed finite volume method in the context introduced in 
|15| using the Raviart-Thomas mixed finite element and this 'broken' Pi -nonconforming element. The advantage of this 
mixed finite volume method is that once we solve the symmetric positive definite pressure equation(without Lagrangian 
multiplier), the velocity can be computed locally by a simple formula. This procedure avoids solving the saddle point 
problem. Furthermore, we show optimal error estimates of velocity and pressure in our mixed finite volume method. 
Numerical results show optimal orders of error in L 2 -norm and broken P 1 -norm for the pressure, and in .ff (div)-norm 
for the velocity. 
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1. Introduction. There are many physical problems where the underlying partial differential 
equations have an interface. For example, second order elliptic equations with discontinuous coeffi- 
cients are often used to model problems in material sciences and porous media when two or more 
distinct materials or media with different conductivities, densities or permeability are involved. The 
solution of these interface problems must satisfy interface jump conditions due to conservation laws. If 
the interface is smooth enough, then the solution of the interface problem is also smooth in individual 
regions where the coefficient is smooth, but due to the jump of the coefficient across the interface, the 
global regularity is usually low and the solution usually belongs to H 1+a (£l) for some < a < 1 . Be- 
cause of the low global regularity, achieving accuracy is difficult with standard finite element methods, 
unless the elements fit with the interface of general shape. 

Immersed interface method using uniform grids has many advantages over usual fitted grid 
method. Using uniform grid, one does not need to generate a grid. This is quite convenient in 
several aspects. First of all, the structure of stiffness matrix is the same as that of the standard 
finite element method, where many known efficient solvers can be exploited. Second, when a moving 
interface problem is involved one does not need to generate a new grid as time evolves. This saves 
considerable amount of time and storage. 

The first attempt to avoid fitted grid for interface problem was made by LeVeque and Z. Li 
|30j . where they proposed an immersed interface method for finite difference method where the jump 
condition was properly incorporated in the scheme. Cartesian grids is most natural in this case. 
They subsequently applied the same idea to other interface problems such as the Stokes flow problem, 
one-dimensional moving interface problem and Hele-Shaw flow, etc. [26j [31] [331 [33] The resulting 
linear systems from these methods are non-symmetric and indefinite even when the original problem 
is self-adjoint and uniformly elliptic. Although these methods were demonstrated to be very effective, 
convergence analysis of related finite difference methods are extremely difficult and are still open. 

For finite element methods, T. Lin et. al. [351 (Ml GUI [35] recently studied an immersed interface 
finite element method using uniform grids and they proved the approximation property of the finite 
element space of their scheme. Their numerical examples demonstrated optimal orders of the error. 
Other related works in this direction can be found in [TOj, ESI Ell [3H [39l S3] and references therein. 
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On the other hand, Pi-nonconforming finite element method introduced in [20] for solving Stokes 
equation is being widely used in solving elliptic equations and shown to be quite effective [15l [TBI l20l 
127) . Especially, it is extremely useful in solving mixed finite element method by hybridization [TJ [2] 
or finite volume formulation [15j HHJ \T7\ [19] . 

The mixed finite element method based on the dual formulation is well-known [51 El EZl El (HI 
|2"T1 |2"31 |2~41 W2\ . The motivation of mixed method is to obtain an accurate approximation of the flow 
variable and has been widely used in the study of flow in porous media such as petroleum engineering, 
underground water flow, and electrodynamics, etc. (e.g., p~2j [22l [41] ) But this scheme leads to a saddle 
point problem for which many well known fast iterative methods fails. To overcome this difficulty, 
mixed hybrid methods have been introduced [2] [8] 140] , where the problem reduces to a symmetric 
positive definite system in Lagrange multiplier only. The flow and pressure variables are obtained via 
some post -processing. 

Recently, there has been some development of the mixed finite element method in another direc- 
tion: A mixed finite volume method was proposed in [19] and extended in |15[ I16j. In this method, 
one use Raviart-Thomas space and Pi -nonconforming space as trial spaces for velocity and pressure, 
and integrates the mixed system of equation on each volume. Then one can eliminate velocity variable 
and obtain the equation of pressure variable only (in terms of Pi-nonconforming FEM) directly from 
the formulation without using Lagrange multiplier. The resulting linear system is again symmetric 
positive definite, and velocity can be recovered from pressure locally in a simple manner. 

The purpose of this paper is two-folded. First, we propose a finite element method on a uniform 
triangular grid using 'broken' Pi-functions having degrees of freedom on edges. This is a Galerkin 
type Pi -nonconforming finite clement method with the basis functions having the average on edges 
as degrees of freedom, broken along the interface to match the flux condition. Then we show optimal 
error estimates in H , L 2 -norms. Here, we emphasize that the meaning of 'nonconforming' is different 
from the context of Li et al. [37l [35] where the basis function has degrees of freedom at vertices, 
discontinuous along edges of interface elements. Meanwhile, the basis function here are Crouzeix- 
Raviart type [20] . Hence it is discontinuous along all edges intrinsically. Furthermore, since we use the 
average of linear function(possibly broken) along edges as degrees of freedom, the overhead of dealing 
with nonconformity in the proof of error estimate is significantly reduced. (See section 3) 

Next, we propose a mixed finite volume method using Raviart-Thomas space and the immersed 
interface finite element introduced above. This is similar to the scheme studied in [15] . but the usual 
nonconforming basis function is replaced by a 'broken' one on the interface element. We provide an 
optimal error analysis of pressure and velocity. 

The rest of the paper is organized as follows. In the next section, we will describe the model 
problem and some preliminaries. We construct an immersed interface Pi-nonconforming space with 
average degrees of freedom which preserves flux continuity weakly along the interface, and prove an 
interpolation error estimate. In Sections [3] and [¥[ we propose an immersed interface finite element 
scheme and prove H 1 and L 2 -error estimates. In Section [5] we propose a mixed finite volume method 
using Raviart-Thomas mixed finite element and our Pi -nonconforming immersed interface finite ele- 
ment method, where the problem reduces to symmetric positive definite system in pressure variables. 
The velocity can be computed locally after pressure computation. Finally, in Section[6l some numerical 
results are presented which indicate optimal orders convergence of our methods. 

2. Preliminaries. Let fi be a convex polygonal domain in K 2 which is separated into two sub- 
domains £1 + and f2~ by a C 2 -interface T = <9f2~ C ft, with £1+ = f2 \ as in Fig. 12.11 We consider 
the following elliptic interface problem 



-div(/3Vp) = / in ft \ T, 
p = on <9f2 



(2.1) 
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Fig. 2.1. A sketch of the domain f! for the interface problem 

with the jump conditions on the interface 

[p]=0, [j9^]=0 across T, 
an 



(2.2) 



where / £ L 2 (S1) and p 6 TJg (CI). We assume that the coefficient f3 is positive and piecewise constant, 
that is, (3{x) = (3- for x en'; 0(x) = (3+ for x £ fl+ . 

We take as usual the weak formulation of the interface problem: Find p £ Hq(CI) such that 

/ P\7p-Vqdx= I fqdx, VqeH^(Sl). (2.3) 
Jq Jn 

Now we introduce the space 

H 2 (fl) := {p £H l (fl) :pe H 2 {fl s ), s = +, -} 

equipped with the norm 

lbll| 2(0) := \\pf HH a) + lbll^ ( n+) + \\pf H *(a-), G ^ 2 (^) ; 

where H m (il s ) — W™(Cl s ) is the usual Sobolev space of order m. By Sobolev embedding theorem, 
for any p £ H 2 (Cl), we have p £ W}(ft), Vs > 2. Then we have the following regularity theorem for 
the weak solution p of the variational problem (|2.3[) ; see [4] and [28] . 

Theorem 2.1. XTie variational problem \2.S\) has a unique solution p £ H 2 (Q) which satisfies 
for some constant C > 



|b||j}2 ( Cj) < C||/|| L 2(fi). 



(2.4) 



We now describe an immersed interface finite element method with piecewise Pi -nonconforming 
functions. 

For the simplicity of presentation, we assume that fl is a rectangular domain. First we consider 
uniform rectangular partitions of mesh size h. Then we obtain triangular partitions 7^ by cutting the 
elements along diagonals. Thus we allow the interface T to cut through the elements. We assume the 
following situation: The interface 

• meets the edges of an interface element at no more than two points. 

• meets each edge at most once except possibly it passes through two vertices. 
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These assumptions are reasonable if we choose h sufficiently small. 

We call an element T G Tj, an interface element if the interface V passes through the interior of 
T, otherwise we call T a non- interface element. (If one of the edges is part of the interface, then the 
element is a non-interface element.) Let £h be a collection of all edges of Th- 

Let DE be the line segment connecting the intersections of the interface and the edges of a triangle 
T. This line segment divides T into two parts T+ and T~ with T = T + U T~ uM There is a small 
region in T such that T r = T - (0+ n T+) - n T~)(see figured). Since JJE can be considered 
as an approximation of the C 2 -curve L n T, the interface is perturbed by a 0(h 2 ) term. From [4lll3|. 
one can see for the interpolation polynomial defined below, such a perturbation will only affect the 
interpolation error to the order of h 2 . 

As usual, we want to construct local basis functions on each element T of the partition Th . For a 
non- interface element T £ Th, we simply use the standard linear shape functions on T having degrees 
of freedom at the mid-points of the edges, and use Sh{T) to denote the linear spaces spanned by the 
three nodal basis functions on T: Let nn, i = 1, 2, 3 be the midpoints of edges of T. Then 

Sh(T) — span{ fa : fa is linear on T and fa(rrij) = &y, i,j = 1, 2, 3 } (2-5) 

Alternatively, we can use average values along edges ej of T as degrees of freedom, i.e., fa can defined 

h y ]k\ Lj fa ds = S v ' *' 3 = 1 ' 2 ' 3 - 

For this space, we have the following well-known approximation property [181 120] : 



lb - hp\\L*(r) + h\\p - Ihp\\m{T) < C^WpWh^t), (2.6) 

where It ■ H 2 (T) — » Sh{T) is the interpolation operator. Finally, we use 5^(0) to denote the space 
of the standard piecewise Pi -nonconforming space with vanishing boundary nodal values. 

2.1. Local basis functions on an interface element. We now consider a typical reference 
interface element T whose geometric configuration is given in Fig. 12.21 in which the curve between 
points D and E is part of the interface. Let e,, i = 1, 2, 3 be the edges of T. For e iJ 1 (T), let 4> ei 
denote the average of <j) along e,, i.e., (p ei := jj-j J e tfids. 

We construct a piecewise linear function of the form 

rh(Y\ - / =ao + boX + c y, X = (x,y)eT+, , 

n *>~ \ 4>-{X) = a 1 +b 1 x + c 1 y, X=(x,y)eT-, 

satisfying 

fa z =Vi, i = 1,2,3, (2.8) 
0+(U) - r(D), ^ + {E) = tf>-{E), (3+^- = ir-j^—, (2.9) 

where Vi, % = 1,2, 3 are given values and n-^ is the unit normal vector on the line segment DE. This 
is a piecewise linear function on T that satisfies the homogeneous jump conditions along DE. 

Suppose that a typical reference interface element T has vertices at A(Q, 0), B(l,0), C(0, 1). We 
assume that the interface meets with the edges at D(xq,0) and E(0,y o ) where < £0,2/0 < 1- Then 
the unit normal vector to the interface is n-^ = (j/o, x o)l V x o + Vo- 

Theorem 2.2. Given an reference interface triangle, the piecewise linear function fax, y) defined 
by \2.7\)-$2lfy is uniquely determined by three conditions 



fa^Vi, i = 1,2,3. 




Proof. Let X = (x, y) T G T . Since (j> + and tj> are linear functions, we have 

4>{X)-- 

The condition (|2.8p gives the following three equations: 



(j) + {X) = a Q + b Q x + c Q y, XeT+, 
<f>~(X) = a x +b lX + c iy , XeT-. 



ao + 7^0 + ic = Vi 



and 



>ds = 



yo 



ds 



AE 



EC 



= («i + y c i)2/o + («o 



yo + 1 



ds 



co)(i - yo) = v 2 . 



where we used mid-point quadrature on AE and EC. Similarly, we have 

4> e3 = (Ol + Y bl ^ X ° + (° + ^2 _ ^ ) = ^ 

From the continuity condition at Z) and -E, we have 

a + b x = oi + bix , 
ao + co2/o =01 + cit/o 
and the flux continuity condition along DE gives 

(&Q)Co) • (yo,x ) = p(h,ci) ■ (yo,x ), 



(2.10) 



(2.11) 



(2.12) 



(2.13) 



(2.14) 
(2.15) 



(2.16) 



where p = [3~ / f3 + and we have used that the normal direction of the line segment DE is (yo,Xo). 

Then the coefficient matrix of the above linear system for the unknowns ao,bo, cq and oi, b\, c\ in 
this order is 



A 



f 1 


1 

2 


i 








\ 


i - yo 





|(i -vl) 


yo 





bl 


1 - x 


\{l-xl) 





Xq 


2 X 





-1 


-Xq 





1 


x 





-1 





-yo 


1 





yo 


V o 


-yo 


-X 





pyo 


P^o / 



(2.17) 
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Tedious calculation(see appendix) shows that the determinant of the matrix is 

det(A) = i(a^ + yl){p{x yo - 1) - x y } < 0. (2.18) 

Thus the coefficients of (|2.10l) arc uniquely determined. □ 

Remark 2.1. If <f> ei ,(j) e2 and <f> e3 have the same value, then the piecewise linear function cj> sat- 
isfying 12.8\) - Un%) reduces to a constant by uniqueness. 

Now we can construct nodal basis functions on an interface element T in general position through 
affine mapping. We let Sh{T) to denote the three-dimensional linear space spanned by these shape 
functions. We note that Sh{T) is a subspace of H (T). Finally, we define the immersed interface 
finite element space Sh(&) as the collection of functions such that 

t £ Sh(T), if T is a noninterface element, 
t £ Su{T), if r is an interface element, 

J ds = J (j)\T 2 ds 7 if Ti, T% are adjacent elements and e is a common edge of T\ and Ti, 
ds = 0, if e is part of the boundary dft. 

Although for functions in Sh(T) the flux jump condition is enforced on line segments, they actu- 
ally satisfy a weak flux jump condition along the interface. This is stated in the following lemma [36], 
whose proof is a simple application of the divergence theorem. 

Lemma 2.3. For an interface triangle T, every function (p G Sh(T) satisfies the flux jump 
condition on T (~)T in the following weak sense: 

[ (/rV(?r -/3+V0+) -n r ds = 0. (2.19) 
JrnT 

Proof. Let <f> be any function in Sh{T). By the divergence theorem, we have 

/ (/rV<T - f3 + \7(j) + )-n r ds+ f(p-V<tr ~[3 + V(f) + ) ■ n-^ds 
JrnT Jde 

div(/T - f3 + V4> + ) dx = 0. 



By the flux continuity of <j) on DE, 



[_ 08" - /3+ V0+) • n^ds = 0, 
Jde 



which completes the proof. □ 



2.2. Approximation property of nonconforming immersed interface space Sh(T). In 
this subsection, we would like to study the approximation property of Sh(T) by defining an interpo- 
lation operator. The difficulty lies in the fact that Sh(T) does not belong to H 2 (T), the restriction 
of H 2 (fl) on T, where H 2 (T) = iJ x (T) n H 2 (T n 0+) n H 2 (T n Q-)(see Fig. H3J). To overcome the 
difficulty, we introduce a bigger space which contains both of these spaces. 



Pi-Nonconforming Finite Element Method For Interface Problems 



7 



For a given interface element T, we consider a function space X(T) such that every p G X(T) 
satisfies 

( pe h\t) n h 2 (t + n n+) n f 2 (t^ n fr) n iJ 2 (r+) n h 2 (t-), 

I / (/T Vp~ - /3+Vp+) • n r ds = 0, 
I JrnT 

where = T r D s , s = +, — . For any p G X(T), we define the following norms. 

Hx(T) = bb.r-nfi- + bl2,T+nn+ + H^t- bl^r^' 

\\p\\x(T) = \\P\\1,T + \p\x(T)> 
3 

|M|2,T = \P\X(T) + ^\PeX 
i=l 

where p ei , i = 1, 2, 3 are the average on each edge e^. 




(a) S h (T) C H 2 (T+) n H 2 {T~) (b) H 2 (T) C if 2 (T n n ff 2 (T n Q - ) 

Fig. 2.3. Regions of H 2 -regularity 

Remark 2.2. If p £ H 2 {VL), then p\ T G X(T) and |p|x(r) = blg2(r)> w/iere blj^m = 
bla.Tnn- + bl2,Tnn+" 

Lemma 2.4. ||| • |||2,t is a norm in the space X(T) which is equivalent to || • ||x(T)- 
Proof. Let p G X(T). If \\p\\ 2 ,T = 0, then |p| x(T) = and \p ei \ = 0, i = 1,2,3. Hence p is 
linear on each of the four regions T + n T~ n Q~, T+ and T°. Since p G i? 1 (T), p is linear on 
each T s , s = +, — . Since p satisfies flux continuity condition, Sh{T). Now we apply Theorem 12.21 to 
conclude p = 0. 

We now show the equivalence of ||| • |||2,t an d || • ||x(T) ( c f- [31 P-77], [HE])- First, note that by 
Sobolev embedding theorem, H 2 (Ti) is compactly embedded in Wg(Ti)fat any s > 2, where Tj C T. 
So we see that X(T) C W}(T) C C°(T). If p € X(T), then p is a continuous function on T and 
|PeJ < C|blU(T), i= 1,2, 3, thus 

|p|a,r < C||p||x(r), (2-20) 

where C is independent of p. 

Now suppose that the converse 

IblUcn < c|p| 2>T , Vp g x(T) 
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fails for any C > 0. Then there exists a sequence {pk} in X(T) with 

\\Pk\\x(T) = h ||k|| 2 ,T < p fc = l,2,-... (2.21) 

Since (T), (s > 2) is compactly imbedded in iJ 1 (T) by Kondrasov theorem [I8j p. 114], there exists 
a subsequence of {pfc} which converges in i? 1 (T). Without loss of generality, we can assume that the 
sequence itself converges. Then {pk} is a Cauchy sequence in 7? 1 (T). Noting that \pk\x{T) — * and 
Ibfc-Pdlx(T) - Ibfe ~P'I!i.t + (bfclx(T) + bz|x(T)) 2 , we see that {p k } is a Cauchy sequence in X(T). 
By completeness, it converges to an element p* £ X(T), and (|2.20p . (|2.2ip gives 



2,T 



< 



Pfclh.T + Ibfelb.T < C||p* -Pk\\x(r) + j. — *• 0- 



But 

lb*llx(r) = 1 and ||b*||| 2 ,T = 0. 
This is a contradiction, since ||b*|||2,T = implies p* = 0. □ 

For any p E X{T), we define I/j.p E Sh(T) using the average of p on each edge by 

(4?)e, =Pe», i = 1,2,3 

and call IhP the interpolant of p in Sh(T). We then define T^j? for p E H 2 (fl) by (I/ip)|t = ^(p|t)- 
Lemma 2.5. Let T 6e an interface element. Then for any p E X(T), we /lave 

lb-4p|U,T<C/i 2 -" l ||p||x(T), m = 0,l, (2.22) 

where h is the mesh size of T . 

Proof. Let T be a reference interface element. Then for any p E X(T) 



4p||| 2i f = |p - 4p| x( f) + ^2 \@ ~ Ih P^ 

i=l 

= \P — IhP\ X (f) = blx(f)' 



where we used the fact that p ei = (Ihfi) ei on each edge and 77 2 -seminorm of the piecewise linear 
function I^p vanishes. Applying the scaling argument for m = 0, 1, we have 

lb - IhP\\m,T < Ch^Wp - I h p\\ m f < Ch l - m \p - I h p\\ 2 f 

< Ch 1 ' m \p\ x(f) <Ch 2 - m \p\ x(T) . 

□ 

By above lemma, Remark 12.21 and (|2.6[) . we obtain the following interpolation estimate. 
Theorem 2.6. For any p E H 2 (fl), there exists a constant C > such that 

lb - IhP\\v<n) + h\\p - hp\\i. h < C/i 2 ||p|| 52(n) , (2.23) 
where \\ ■ \\ 2 h := J2reT h II ' lli.r- 
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3. Immersed interface FEM with 'broken' Pi-nonconforming elements. We are now 

ready to define our immersed interface finite element method based on 'broken' Pi-nonconforming 
element: Find ph £ Sh(&) such that 

a h (p h ,<t>h) = {f,<Ph), W<f> h €S h (Q), (3.1) 

where 

a h (p,cj>)= f PVp-Vct>dx, Vp,<peH h (il), (3.2) 

H h (Q) := H^(Q)+S h (Q). 

Here, Hh(Q) is endowed with the piecewise iJ 1 -norm || • Hi^. Note that if discrete Poincare inequality 
holds, then noting that the bilinear operator a/j(-, •) is bounded and coercive on Sh(Q), the discrete 
problem (|3.1|) has a unique solution ph G Sh (^) ■ 

Lemma 3.1 (Discrete Poincare inequality). There exists a constant C > independent of h such 
that for any (f> £ SVi(^) 

C\\4>\\l, [a) <a h (<f>,<t>). (3.3) 

Proof. Let e be the common edge of two adjacent elements T\ and T%. Note that since J <j>i ds = 
J e <p2ds, where <f>i = </>|t (! i = 1,2, there exists a point io £ e such that 4>i{ x o) — 4>2{xq). Then a 
slight modification of Lemma 2.1 in 14] proves the inequality. □ 

For the energy-norm error estimate of the immersed interface finite element method, we need the 
well-known second Strang Lemma which is valid since a%{-, •) is coercive. 

Lemma 3.2 (Second Strang Lemma). If p e H 2 (Q), ph € S/,(0) are the solutions of &2.3\) and 
\3.1\ ) respectively, then there exists a constant C > such that 

ii n s n) ; t II ii I ah{p,4>h) ~ UvM I [ , Q 
\\p-Ph\\i.h < C< lnf \\p-qh\\i,h+ sup [7— -jj }. (3.4) 



We shall need the following estimate; see Lemma 3 in [20] . 

Lemma 3.3. Let e be an edge of T. Then there exists a constant C > such that for all 



(v — v e ) ds 



< Ch\<f>\i, T \v\i, T) 



where v e '■= -rjy f e v ds 



Remark 3.1. This lemma also holds when <p belongs to H (T s ), s = =t with |<^|i,t understood as 
sum of piecewise norm |0|ix±- 

Theorem 3.4. Let p e H 2 (fl), p h e Sh(Q) be the solutions of S2.3\) and \3.1\) respectively. Then 
there exists a constant C > such that 

\\p-Ph\\i,h < Ch\\p\\ R2{ny (3.5) 
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Proof. We use the second Strang Lemma. The first term is nothing but an approximation error. 
By Theorem 12.61 we have 

hif \\p-qh\\i,h < Ch\\p\\g 2 , n y (3.6) 

9(i6S h (n) 

For the consistency error, we have from the definition of a^-, •) and Green's formula 

a h (p,(j) h ) - (f,(j)h) = Y / PVp-\7(j) h dx - / f4> h dx 
TeT h Jt ^° 



dT 

(hi 

TeT h 

= E <^^> 9 T, (3.7) 

TeT h 

where 4>h E Sh{Q) and n is a unit outward normal vector on each dT. Since /?§f belongs to ff x (Tn 
nif^Tnfi - ) and (f>h E Sh(£l) has well-defined average value on the interior edges, and vanishing 
average on the boundary, we have by Lemma 13.31 and remark 13.11 



E<^*>-EE<^-(^)..*>. 

TeT h TeT h eCST 

^ E ^I/3^Ii,t|^|i,t 
TeT h n 

<Ch\\p\\^ m \\ct> h \\i, h . (3.8) 

This completes the proof. □ 

4. L 2 -error estimate. We now apply the duality argument to obtain L 2 -norm estimate of the 
error. Let us consider an auxiliary problem: Given g E L 2 (Q), find (p £ H 2 (fl) such that 

- div(/3V<y9) = g mO\r, (4.1) 
ip = on dfl 

with jump conditions [u] = 0, [/3§^ ] = across F. Then we have 

IMIj? 2( n) <C\\g\\v [a) . (4.2) 
Let ifh E S/i(f2) be the solution of the corresponding variational problem 

ahivhiVh) = {v h ,g), Vv h eS h (n). (4.3) 

Then 

(p-Ph,9) = Y j P^{P~Ph) -^fdx - Y j (p-Ph)P^ds 
TeT h T TeT h dT n 

= a h (p ~Ph,<P - fh) + a h (p -p h , f h ) - Y / iP-Ph)P-~-ds 

TeT h , ' dT 71 

= a h {jp-p h ,ip-ip h )+ Y j P^^hds-Y ( (p-Ph)/3^r-ds 
TeT h dT n TeT h dT U 

=: I + II - III. 
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By continuity of a/j(-, ■) and /Z^-error estimate of ip — <ph, 



\I\ < C\\p - p h \\ h h\\<p - <Ph\\i,h < Ch\\p - ph\\i,h\\(p\\H2( n) 

< Ch 2 \\p\\ 52(Q) \\(p\\ S2(n) . 

Applying the analysis for the consistency error of i/^-error estimate ()3.8j) . we get 



\II\ 



Y] [ P^-Vhds = V / /3^((p h -tp)ds 



TeT h 



< Ch\\p\\ S2{n) \\<p h - (fWu, < C^WpW^^WvWhuq) 



and 



\M\ = 



52 / (p-Ph)P^rds 

hi. Jot dn 



T&T h 
2 



< Ch\\p - p h \\i,h\\<p\\ff2 



(£2) 



< Ch ||p||52 (n) ||v||g 3(Q) . 
Since ||v|[g2(n) < C||.9lU 2 (o), we see that 



\\p-Ph\\L'W= sup -j-J^- < C^IHIjp 



g eL 2 (n) \\g\\L 2 (n) 
Thus we obtain the following L 2 -error estimate. 



(ay 



(4.4) 



Theorem 4.1. Letp G H 2 (tt), p h G <S/i(fi) be the solutions of \2.3\) and (3.1\) respectively. Then 
there exists a constant C > such that 



\\P-Phh»(n) < Ch 2 \\p\\ M2m . 



(4.5) 



5. Mixed finite volume method based on IIFEM. In this section, we propose a new mixed 
finite volume method based on the 'broken' Pi-nonconforming interface finite element method intro- 
duced in the previous section. Our method is similar to the mixed finite volume method studied 
in [TBI H"6l [19], but the usual nonconforming finite element space is replaced by our 'broken' P\- 
nonconforming space. 

Let us write the problem (|2.3p in a mixed form by introducing the vector variable u = — /3Vp as 



u + (3Vp = in fi, 
divu = / in 17, 
p = on dft. 



(5.1) 



The mixed finite element method based on this dual formulation is well-known [7J 32] • The idea 
of the mixed method is to find a direct approximation of the flow variable u. For that purpose, we 
introduce V = H(div, f2) = {v G L 2 (il) : divv G i 2 (fi)}, and use the local RTq space to approximate 
the flow variable which is given by Vft(T) = {v : v = (a + ex. b + cy), a,b,c G M} for any triangle 
element T. The global space V?, is defined as 



Vft, = {v : v\t G Vfe(T); v • n is continuous along interior edges}. 



(5.2) 
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This method gives a good approximation of the flow variable. However, it leads to a saddle point 
problem, that is, one obtains an indefinite matrix system when (|5.ip is discretized. As mentioned 
earlier, a popular way to avoid this indefinite system is to use Lagrange multipliers [2]. Another 
possibility is to form a mixed finite volume method as in pj2 UH [19] . 

To define a mixed finite volume method for an interface problem, we use the well-known RTq space 
V/j for velocity and 'broken' Pi -nonconforming immersed interface space Sh for pressure variable. Note 
that every v 6 has continuous normal components across the edges of Th, which are constant. 

We consider the following scheme: Find (uh,Ph) G V/i x Sh which satisfies on each element T G Th 

f {u h + pVp h ) ■ V<t> = 0, V0 e S h (T), 

r r (5-3) 

/ divuh = / /. 
Jt jt 

Note that since divu^ is constant, divuh\r = f\r '■= wr J T /, where |T| denotes the area of T. When 

the interface is not present, Sh(T) = Sh(T) and this scheme coincides with the one in [15l [T9] . Since 
the numbers of unknowns and equations do not change, our scheme is a square linear system and has 
a unique solution. We refer to [15] for details. 

Now since • n is constant on the edge and (f) G Sh has common average values on interior edges 
and vanishing boundary nodal values, we obtain 

]T f u h • V0 = f ( u * • - / divu ^ = - / 70, (5.4) 

where / G L 2 (fl) is a simple function having value f\r for each T. From (|5.3p . it immediately follows 
that 

Y, f /3V P/i .V0= f /0, Vcf>eS h (n). (5.5) 

This is the interface finite element method introduced in the previous section, except that on the 
right-hand side / is replaced by /. 

The velocity can be computed directly from the solution ph of (|5.5|) as follows. Let T be an 
any element of Th with the edges e^, i = 1, 2, 3 and let <pi G Sh(T) be the 'broken' Pi-nonconforming 
basis function associated with the edge ei. Then the flux through the edge is given by 

\ei\(u h -n)| ei = / (u h -n)4>i = / div(u /l 4 ) 

JdT JT 

= / (divu/,^ + u ?l • V^i), 

JT 

where 4>i G Sh(T) is a basis function on T. Then it follows by (|5.3p that 

|ei|(u fc • n)| ej = / - / pV Ph ■ V^. (5.6) 

JT JT 

Thus in order to compute the fluxes through the edges of an element T, we only need to compute the 
local residual of the solution ph on each T. 

The error estimate of U/j would follow that of ph- In fact, we can relate the estimate ||u — U/, || 
with ||p — PfcHi,/^ First, we show the following local formula. 



Lemma 5.1. Let u/j, ph be the solutions of \5.3\) , then 

7, 



u h (x) = -/3Vp h + |(x-x B ), VxgT, (5.7) 
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where fJVph denotes the average of (3\7ph on T and xg is the center of T . 
Proof. Expanding about x^, the barycenter of T, we have 

Ufc(x) = u h (x B ) + T>u h (x B )(x-x B ), x e T, (5.8) 

where Puj, is the Jacobian matrix of u^. Let u/j = (a + cx, b + cy) 6 V^(T), then we have 

Vu h (x B )(x - x B ) = c(x - x B ) = dl ™ h (x - x B ) = -(x - x B ), 

where we used the relation divu/j = /. On the other hand, applying cf) — (x, 0) T and (0, y) T in (|5.3[) . 
we see 

-WPh = T^Ju h = u h (x B ). (5.9) 
Substituting these into (|5.8[) . we obtain formula (|5.7j) . □ 

Remark 5.1. Our formula is different from the one in U51 \16\j where they have 

u h (x)=-/3V Ph + £(x-x B ), VxeT. (5.10) 

The ph in our scheme is broken along a line segment contained in an interface element, Hence we 
have taken the average of (3Vph . 

By the above lemma, we have 

7 

u(x)\ T - u h (x) = -f3Vp + f]Vp h --(x-x B ). 

So 

||u-Ufc|| 0j T < ||/3Vp-^V^||o,r + C|7|||x-x B || ,r. 
Since (3 is piecewise constant and u = — /3Vp, we have 

||/3Vp - Wp1\\o,t < H/SVp - Wp\\o,t + \\(3Vp - l3V Ph \\ ,T 

< Ch\\u\\i tT + Ch\(3V P - (3W Ph \ 

< Ch\\u\\x t T + C\p - Ph\l,T, 

provided u G H 1 (T). Hence 

||U - U h || ,T < C{h\\u\\ hT + \p- p h \ hT + h 2 \J\} 

< C{h\\u\\ hT + \p - Ph\l,T + h\\f\\ ,T}, 

where we used |/| < h^ 1 ||/||o,t and ||x — x B ||o,T < Ch 2 . 
Summing over every T £ T^, we have 

||u-Ufc|| L 2 (n) < C{/i||u|| H i(n) + \\v-Ph\\i,h + h\\f\\ L 2 (n) }. (5.11) 

Since divu = / and divu^ = /, we can easily obtain the estimate for ||divu — divu/ l || i 2(r2). This is 
summarized in the following. 

Theorem 5.2. Let u; l; p^ be the solutions of i5.3\) , then there exists a constant C > such that 
u h \\ L 2 {n) + ||divu-divu h || L 2 (n) < C7j{||u]] H i(n) + WpWh^u) + \\f\U,h}, (5-12) 
provided u 6 H 1 (il). 
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6. Numerical examples. In this section, we report numerical results for the schemes introduced 
previously. For a numerical test, we solve problem (|2 . 1 1) with the rectangular domain f2 = [— 1, 1] x 
[— 1, 1] partitioned into unform right triangles having step size h. We take a circle with radius ro = 0.5 
as an interface, and the exact solution is chosen as 




J_. 

/3+- 



in Q 



n+. 



We note that this example is taken from Z. Li [36]. We present errors in L 2 ,iJ 1 -norm for the 
pressure p, while in i?(div)-norm for the velocity u. Here the order of convergence is determined 
by the least squares fit to the data. In Table l6~T1 and HT27 first two columns), we report the results 
of the 'broken' Pi -nonconforming immersed interface scheme introduced in Section [3j where we used 
conjugate gradient method(CG) to solve the resulting discrete system. It shows optimal order of 
convergence for L 2 -norm and if -norm: 



\\p-Ph\\o~0(h 2 ), \\p- Ph \\i, h *0(h). 



(6.1) 



We also present some result for mixed finite volume method introduced in Section [5] Again, this 
shows optimal order of convergence for the flow variable which is consistent with Theorem 15.21 (cf. 
last columns of Table 16.11 and 16. 2p : 



|u - u h \\ L 2 {n) + ||div(u - u h \\ L 2 {n) W 0(h). 



(6.2) 



This is in good agreement with some fitted grid computation; when the jump of the coefficient is large, 
one usually have O(h) order accuracy, see [9] problem 1, p. 310, for example. 



1/h 


\\p-Ph\\o 


order 


\\p-ph\\l,h 


order 


||u - Uh\\o 


order 


||div(u - u h )\\ 


order 


8 


9.576e-3 




1.208c-l 




2. 945c- 1 




1.053c+0 




16 


2.666e-3 


1.845 


6.744e-2 


0.841 


1.702e-l 


0.791 


5.292e-l 


0.993 


32 


6.488e-4 


2.039 


3.341e-2 


1.013 


8.906e-2 


0.934 


2.650e-l 


0.998 


64 


1.400e-4 


2.212 


1.657e-2 


1.012 


4.290e-2 


1.054 


1.326e-l 


0.999 


128 


3.716e-5 


1.914 


8.242e-3 


1.008 


2.015e-2 


1.090 


6.629e-2 


1.000 


256 


8.973e-6 


2.050 


4.117e-3 


1.001 


9.865e-3 


1.030 


3.315e-2 


1.000 


Order 




2.029 




0.985 




0.994 




0.998 



Table 6.1 

Nonconforming immersed interface FEM: j3~ = 1, (3 + = 1000 
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1/h 


\\p-Ph\\o 


order 


\\p-Ph\\i,h 


order 


ll u - Ufello 


order 


div(u - u h )\\ 


order 


8 


1.447e-2 




6.575e-l 




3.361c-l 




1.053c+0 




16 


3.497e-3 


2.049 


3.312e-l 


0.989 


1.657c-l 


1.020 


5.292e-l 


0.993 


32 


8.826e-4 


1.986 


1.661c-l 


0.996 


8.165e-2 


1.021 


2.650e-l 


0.998 


64 


2.210e-4 


1.998 


8.311e-2 


0.999 


4.075e-2 


1.003 


1.326e-l 


0.999 


128 


5.507e-5 


2.005 


4.157e-2 


0.999 


1.959e-2 


1.057 


6.629e-2 


1.000 


256 


1.370e-5 


2.007 


2.079e-2 


1.000 


9.658e-3 


1.020 


3.315e-2 


1.000 


Order 




2.005 




0.997 




1.024 




0.998 



Table 6.2 

Nonconforming immersed interface FEM: f3~ = 1000, /3 + = 1 



Appendix: Computation of determinant of A. 

By adding the last three columns to first the three, we obtain 



1 

2 


W-v 2 ) 













1 


1 

2 


1 

I 

2 














y 





h 2 




1 





y 





\v 2 


~x 2 ) 





X 


l r 2 
2 X 







1 


1 

2 





X 


1-2 
2 





—x 





1 


X 
















1 


X 








-y 


1 





y 













1 





y 


-y 


—x 





py 


px 







(p-l)y 


(p - l)x 





PV 


px 



Now gaussian elimination gives 



1 


1 


1 

2 













1 


1 


1 

2 














\ 

2 





y 





b 2 







\ 

2 





V 





b 2 








1 

2 


X 


1™2 
2 X 













1 

2 


X 


1J 

2 th 














1 


X 
















1 


x 














1 





y 













1 





y 







(p - l)x 





py 


px 










(p - l)x 


2(p-i)y 2 


py 


c 



where C — px + (p — l)y 3 . Continuing 



1 


1 


1 

2 














\ 

2 





V 





h 2 








1 

2 


X 


1-2 
2 














1 


X 














1 





y 











.4 


B 


c 



where 

A = 2(p-l)(x 2 +y 2 ), B = py+(p-l)x 3 . 

Hence the determinant is 

^(xyA -yB- Cx) = -{2{p - l)xy{x 2 + y 2 ) - y{py + {p- l)x 3 ) - x{px + {p- l)y 3 )} 

= -^(x 2 + y 2 ){p( x v - !) - x v}- 



This verifies (|2~TS1) . 
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